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ABSTRACT 

^  All  cost  estimates  are  characterized  by  some  uncer¬ 
tainty.  A  device  helpful  in  communicating  this  uncertain¬ 
ty  to  the  decision  maker  is  a  subjective  probability 
distribution  of  the  system  cost.  A  technique — termed  the 
Subjective  Probability  Estimation  Technique  (SPET) — is  de¬ 
scribed  and  a  computer  program  is  presented  to  facilitate 
its  use.  This  technique  permits  the  analyst  to  represent 
his  notions  about  cost  uncertainty  with  the  beta  or  other 
statistical  distributions. 
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INTRODUCTION 


The  cost  analyst  is  faced  with  many  uncertainties  as 
he  attempts  to  estimate  the  costs  associated  with  a  new, 
undeveloped  system.  He  may  wonder: 

(1)  Will  the  physical  characteristics  of  the  system  re¬ 
main  unchanged  as  the  development  process  proceeds? 

(2)  Will  there  be  any  unforeseen  problems  in  the  develop¬ 
ment  process? 


(3)  Will  the  economic  state  of  the  firms  or  industry  re¬ 
sponsible  for  system  development  and  production  con¬ 
tinue  to  change  as  forecasted? 

(4)  Is  the  quality  of  the  historical  cost  data  suffi¬ 
ciently  high  to  inspire  confidence  in  the  estimates 
made  with  it? 

(5)  Have  the  cost-estimating  relationships  been  properly 
specified? 

To  the  extent  the  analyst  is  unable  to  obtain  complete  an¬ 
swers  to  these  and  similar  questions,  his  cost  estimates 
will  be  enshrov ieci  wiub  unceirtein'ty •  Since  i t  is  impose i*- 
ble  to  obtain  definitive  answers  to  all  these  questions, 
his  cost  estimates  will  always  be  characterized  by  some 
uncertainty. 

The  analyst  can  treat  this  uncertainty  in  one  of  sev¬ 
eral  ways.  He  can  choose  to  ignore  it  and  simply  report 
to  the  decision  maker  the  estimate  which  represents  the 
"most  likely"  or  "best"  estimate  of  cost,  as  in  the  case 
of  this  hypothetical  guided  missile  system: 

FIGURE  1 

"BEST"  UNIT  COST  ESTIMATE  OF  H  HYPOTHETICAL 
GUIDED  MISSILE  SYSTEM 
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This  approach,  however,  belies  the  existence  of  a  range  of 
possible  costs;  when  the  system  is  finally  acquired,  any 


one  of  the  innumerable  costs  within  this  range  may  have 
been  realized.  Such  an  oversimplification  may  mislead  the 
decision  maker  by  causing  him  to  place  excessive  confi¬ 
dence  in  the  best  cost  estimate.  An  illustration  of  this 
potential  pitfall  is  provided  in  the  following  example: 

Suppose  two  alternative  systems  that  do  the  same 
task  are  being  compared.  Suppose,  too,  that  on  bal¬ 
ance,  the  differences  in  effectiveness,  performance, 
growth  potential,  maintainability  and  similar  consid¬ 
erations  between  the  two  systems  are  small,  so  that 
the  choice  is  primarily  one  of  cost.  Suppose  one 
system  is  estimated  to  cost  $1.25  million  and  the 
other  $1.50  million.  Without  an  indication  of  the 
possible  high  and  low  values,  the  $1.25  million 
alternative  would  be  the  logical  choice.  But  the 
choice  becomes  more  difficult  if,  as  shown  in  the 
accompanying  tabulation,  the  $1.25  million  cost  is 
qualified  with  an  estimate  of  a  possible  high  value 
of  $2.00  million,  whereas  for  the  other  alternative 
the  limits  are  estimated  to  be  tighter,  and  the  pos- 


sible  high  value  is 

only  $1.60 

million.  The 

case  for 

Cost  Estimate 

System 

(Millions  of  Dollars) 

-•*  - 

-  -  -  Most 

— 

Lowest 

Likely 

Highest 

First  Alternative 

1.00 

1.25 

2.00 

Second  Alternative 

1.40 

1.50 

1.60 

the  alternative  with  the  "most  likely"  cost  of  $1.25 
million  is  now  more  dubious,  because  its  uncertainty 
spread  is  greater,  extending  on  the  high  end  to 
greater  costs  than  the  other  alternative. 3/ 


The  analyst  can  avoid  the  problem  associated  with  a 
single  best  cost  estimate  by  supplying  the  decision  maker 
with  estimates  of  the  lowest  and  highest  possible  costs  in 
addition  to  the  best  estimates: 


3J  W.  Suthcxland,  Fundamental*  ot  Coit  UnceAtatnty  Analy- 
ilA  (McLean,  VAl  Ra4e.aA.ch  Knaly4lh  Cotip.,  RAC-Ck-4 
777;),  pp  3-4. 
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FIGURE  2 . 

"BEST"  AND  RANGE  ESTIMATE  OF"  THE  UNIT  COST 
OF  A  HYPOTHETICAL  GUIDED  MISSILE  SYSTEM 


’  ?60K  $100K 

This  approach  has  merit  in  that  it  reflects  the  range  of 
the  cost  uncertainty.  However,  it  gives  little  Informa¬ 
tion  about  the  nature  of  the  uncertainty,  e.g.,  whether 
all  the  values  in  the  range  are  almost  equally  likely  to 
occur,  or  whether  the  values  closer  to  the  best  estimate 
are  much  more  likely  to  occur  than  those  near  the  extremi¬ 
ties.  ‘Such  knowledge  could  be  valuable  to  the  decision 
maker. 
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One  device  helpful  in  communicating  knowledge  of  both 
the  range  and  nature  of  the  uncertainty  is  the  probability 
distribution  of  the  total  system  cost.  To  illustrate, 
consider  a  probability  distribution  of  the  cost  of  the  hy¬ 
pothetical  guided  missile  system: 


FIGURE  3 


PROBABILITY  DISTRIBUTION  OF  THE  UNIT  COST 
OF  A  HYPOTHETICAL  GUIDED  MISSILE  SYSTEM 


Note  that  the  range  of,  say,  a  95  percent  confidence  in¬ 
terval  [$80K,  $140K]  is  significantly  smaller  than  the 
full  range  [$60K,  $180KJ .  The  knowledge  that  the  analyst 
is  95  percent  confident  that  the  cost  will  occur  in  the 
much  smaller  interval  [$80X,  $140K]  will  permit  the  deci¬ 
sion  maker  to  act  with  a  more  precise  idea  of  the  probable 
cost  of  the  system  than  he  could  otherwise  (providing,  of 
course,  that  he  trusts  the  analyst's  judgement). 

A  popular  source  for  the  probability  distribution  of 
the  cost  of  a  weapon  system  is  the  prediction  interval  ob¬ 
tained  from  cost-estimating  relationships  (CERs)  developed 
by- regression  analysis.  Unfortunately,  this  method  of 
probability  analysis  has  a  serious  limitation:  there  is 
no  provision  for  the  analyst  to  incorporate  into  the  anal- 
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The  fact  that  this  model  predicts  the  impossible  —  nega¬ 
tive  costs  —  again  reveals  the  above-stated  limitation  of' 
this  approach  to  probability  analysis.  What  the  analyst 
needs  is  a  technique  that  will  permit  him  to  retain  the 
"most-likely"  estimate  of  system  cost  and  incorporate  his 
a  priori  beliefs  into  the  prediction  interval. 

This  paper  describes  such  a  technique.  The  authors 
have  entitled  it  the  Subjective  Probability  Estimation 
Technique  (SPET) .  This  technique  is  based  on  the  same 
principles  used  by  Program  Evaluation  Review  Technique 
(PERT)  analysts  years  ago  to  treat  time-estimating  uncer¬ 
tainty.!./  As  its  name  suggests,  SPET  accounts  for  the 
fact  that  the  probability  distribution  of  the  cost  of  a 
new  system  is  by  necessity  subjective  since  repeated  ob¬ 
servations  on  the  cost  of  the  system,  from  which  an  objec¬ 
tive  probability  distribution  could  be  inferred,  cannot  be 
made  (there  is  only  one  observation  on  the  cost  of  a  new 
system  —  the  final  acquisition  cost  of  the  system  —  and 
when  this  observation  is  made  the  need  for  an  estimate 
terminates) .  The  analyst  implements  SPET  in  three  steps 
by: 

(1)  decomposing  the  system  under  examination  into 
several  subsystems  whose  costs  are  additive; 

(2)  selecting  the  subjective  probability  distribu¬ 
tion  that  best  represents  his  knowledge  and  judgement 
about  the  cost  of  a  subsystem; 

(3)  combining  the  subjective  probability  distribu¬ 
tion  of  each  subsystem  cost  into  a  subjective  probability 
distribution  of  total  system  cost. 

The  remainder  of  this  paper  discusses  the  theoretical  and 
practical  aspects  of  these  three  steps. 

SYSTEM  DECOMPOSITION 

When  a  cost  analyst  estimates  the  cost  of  a  complex 
system  he  generally  breaks  the  system  down  into  several 


F.  S.  Hilticn.  and  G.  J.  Liebcrman,  Introduction  to  Op¬ 
eration*  Research  {San  Francitco ,  CJ1  Hoiden-Vay , 

Inc.  1967 ) ,  pp.  T27  229-232.  See  al&o  K.  R.  Mac Crimmon 
and  C.  A.  Rj/avec,  An  Analytical  Study  oj  the  PERT  A*- 


iumption*  { Santa  Monica 


o 


subsystems  and  estimates  the  cost  of  each  subsystem.  The 
breakdown  of  the  system  is  usually  determined  by  the  ana¬ 
lyst's  knowledge  of  the  system  and  the  form  of  his  data 
base.  Using  the  technique  described  in  this  paper,  the 
analyst  will  also  develop  a  subjective  probability  distri¬ 
bution  describing  his  uncertainty  as  to  the  cost  of  these 
subsystems.  Of  course,  if  some  subsystem  cost  is  known 
precisely,  no  uncertainty  is  involved  and  this  cost  is 
treated  as  a  constant. 

SELECTING  THE  SUBJECTIVE  PROBABILITY  DISTRIBUTION 


A  probability  distribution  can  be  selected  to  repre¬ 
sent  any  imaginable  combination  of  factual  knowledge  and 
subjective  notions  an  analyst  might  have  about  the  cost  of 
a  subsystem.  For  example,  suppose  the  analyst  has  a  good 
idea  of  what  the  lowest  and  highest  possible  costs  for  a 
subsystem  could  be,  but  he  feels  that  all  costs  within 
that  range  are  equally  likely.  His  subjective  probability 
distribution  for  the  cost  of  this  subsystem  can  be  repre¬ 
sented  quite  adequately  with  the  uniform  probability  den¬ 
sity  function: 

FIGURE  7 

UNIFORM  PROBABILITY  DENSITY  FUNCTION 
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Suppose  that  instead  of  feeling  that  all  costs  within 
the  range  are  equally  likely  the  analyst  feels  that  a  par¬ 
ticular  cost  within  the  range  is  more  likely  to  be  real¬ 
ized  than  any  other.  He  could  represent  the  cost  with  a 
triangular  distribution: 

FIGURE  8 

TRIANGULAR  PROBABILITY  DISTRIBUTION 
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or  the  beta  distribution: 

FIGURE  9 

BETA  PROBABILITY  DISTRIBUTION 


or  any  other  distribution  with  a  single  maximum  value. 

The  beta  distribution  is  one  of  the  most  popular 
among  cost  analysts  for  representing  subjective  probabil¬ 
ity  distributions.  The  popularity  of  the  beta  is  due  to 
its  several  appealing  characteristics.  One  of  these  char¬ 
acteristics  is  that  its  range  may  be  restricted  to  posi¬ 
tive  values;  costs  similarly  are  positive.  Another  is 
that  the  beta  has  a  finite,  rather  than  infinite  range;  ic 
is  reasonable  to  suppose  that  the  cost  is  bounded  by  fi¬ 
nite  upper  and  lower  bounds.  Finally,  the  beta  distribu¬ 
tion  can  be  expressed  in  an  infinite  variety  of  skewed  and 
symmetric  forms  which  provide  the  analyst  considerable 
choice  when  specifying  the  particular  shane  of  the  distri¬ 
bution.*/ 


Because  of  the  popularity  of  the  beta  distribution, 
the  discussion  in  the  remainder  of  this  section  will  cen¬ 
ter  on  it.  In  the  next  section,  however,  a  computer  pro¬ 
gram  is  described  that  permits  the  analyst  to  use  any  im¬ 
aginable  subjective  probability  distribution  to  represent 
subsystem  cost. 


5/  A noth.a.n  commonly  uAed  dlAtrlbutlon  Ia  the  Welbull. 

Moa£  oi  the.  appealing  propertleA  o b  the  beta  axe  alAO 
bound.  In  the  W elbull  dlAtrlbutlon.  For  exampltA  o{, 
the  uAe  ob  the  Welbull  In  treating  uncertainty  In  coAt 
analyAlA  Aee  V.  F.  Schaebe.r,  et.  al. ,  A  Monte  Carlo 
Simulation  Approach  to  CoAt- Uncertainty  Analy a xa~, 
(McLean,  l/A :  ReAearch  AnalyAlA  Corp.;  RAC-TP-  349 , 
j 9 6 9 )  and  W.  H.  Sutherland,  A  Method  bo*  Combining 
AAtjmmetrlc  Three-Value  PredlctlonA  oj  Time  or  CoAt 
fMc.Lea.it,  17X1  keAearch  AnalyAlA  Corp. ;  RAC-P-65, 

1972). 
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The  usual  expression  for  the  beta  probability  density 
function  (pdf)  is :!/ 

f (x)  =  Cxa(l  -  x)b;  0  <  x  <  1;  a,b  >  0;  [1) 

where  C  =  r  (a  +  b  +  2)/[T(a  +  l)T(b  +  1)] 

»  » 

=  the  inverse  of  the  complete  beta  function 

and  r  ( t )  =  /°°  zt-le-zdzf  t  >  0. 
o 

This  version  of  the  beta  pdf  will  be  called  the  normalized 
beta  pdf,  since  the  range  of  x  is  the  unit  interval. 

A  simple  linear  transformation,  x*  =  L  +  (H  -  L)x, 
where  L  and  H  are  the  lowest  and  highest  values  of  x*,  re¬ 
spectively,  extends  the  range  of  x  in  equation  [1]  to  any 
finite  interval,  yielding  a  generalized  beta  pdf: 

g  (x*)  =  [C/(H  -  L)a+b+1J(x*  -  L)a(H  -  x*)b  [2] 


where  C  is  as  defined 

a ,  b  >  0  > 


ir.  equation  [1]  and  L»  x*  H, 


The  four  parameters  of  the  generalized  beta  pdf  are  a,  b, 
L,  and  H ,J_/  Therefore,  a  unique  pdf  is  defined  for  every 
four-tuple  (a,b,L,H).  The  values  the  analyst  assigns  to 
these  parameters  can  be  obtained  through  certain  estima¬ 
tion  procedures. 

The  analyst  may  estimate  L  and  H  directly  from  his 
and  other  expert  knowledge  of  the  subsystem's  technology, 
contractor  (builder) ,  industry,  etc.  After  analyzing 


bj  H.  J.  Larson,  Introduction  to  Probability  Theory  and 
Statistical  inference  (New York:  John  Wiley  and  Sons , 
Inc.,  1969),  pY  305;  8.  W.  Lindgren ,  Statistical  Theo¬ 
ry  (New  York:  The  MacMillian  Co.,  19  68),  p,  31 3. 

TJ  Note  that  L  and  H  serve  only  to  specify  the  origin  and 
range  o&  x* ,  uihereas  a  and  b  determine  the  shape  oi 
the  pdf  oi  x*. 
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these  data,  he  chooses  L  and  H  so  that  the  cost  of  the 
subsystem  could  never  be  less  than  L  or  greater  than  H. 

Estimates  of  the  parameters  a  and  b,  however,  cannot 
be  obtained  in  a  direct  manner.  One  way  to  estimate  them 
is  to  obtain  two  functionally  independent  equations  in  a 
and  b  and  then  solve  them  for  these  parameters.  In  Appen¬ 
dix  A,  two  such  equations  are  proposed  and  a  computer  pro¬ 
gram  facilitating  their  solution  is  documented.  Potential 
users  of  this  method  for  estimating  a  and  b  are  cautioned; 
some  combinations  of  the  analyst-supplied  inputs  result  in 
distributions  that  cannot  be  represented  by  a  beta  random 
variable.  This  problem  can  be  avoided  by  using  another 
pair  of  equations.  Consider  the  mode  M(x*)  and  the  vari¬ 
ance  V(x*)  of  the  generalized  beta  distribution :£/ 

M (x*)  =  (aH  +  bL) / (a  +  b) ,  L  =  M(x*)  £  H  [3) 

V (x*)  =  [(a  +  1) / (b  +  1)  (H  -  L)2]/[a  +  b  +  2)2 

(a  +  b  +  3)],  0  £  V (x*)  =  (H  -  L) 2/12  [4] 

Using  cost-estimating  relationships  or  other  methods, 

the  analyst  can  estimate  the  mode  or  most  likely  value  of 
the  subsystem  cost.9/  By  evaluating  much  of  the  informa¬ 
tion  he  has  about  the  subsystem  cost  he  can  estimate  its 


£/  See  Footnote.  I  In  Appendix  A  fo> i  the  derivation  of 
thebe  formuZae. 

The  range  of  the  mode  is  obtained  from  its  definition. 
The  range  ojj  the  variance  is  due,  in  part,  to  the  fact 
that  the  Zotver  bound  on  the  variance  of  any  distribu¬ 
tion  is  zero,  and,  in  part,  to  the  fact  that  the  beta 
distribution  converges  to  the  uniform  as  parameters  a 
and  b  approach  zero.  The  upper  bound  on  the  variance 
of  the  generaZized  beta  distribution  is,  therefore, 
the  variance  of  the  [generaZized)  uniform  distribu¬ 
tion,  nameZy 

iMx*)  «  [H  -  L) 2/ J  Z  [5] 

9/  Under  most  circumstances  the  most-ZikeZy  vaZue  is  the 
anaZyst's  point  estimate  of  the  subsystem  cost. 


.St*."-  *»■„. 
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variance. 12/  Having  obtained  estimates  for  these  two  sta¬ 
tistics,  he  can  solve  equations  [3]  and  [4]  simultaneously 
to  determine  unique  values  (estimates)  for  a  and  b. 

The  principal  difficulty  with  the  technique  proposed 
is  in  the  estimation  of  the  variance.  It  is  difficult  for 
the  analyst  to  interpret  his  beliefs  concerning  the  uncer¬ 
tainty  surrounding  a  point  estimate  in  terms  of  the  beta 
variance.  As  an  aid  in  this  process,  the  analyst  may  con¬ 
sider  a  related  variable,  which  is  termed  an  uncertainty 
coefficient  in  this  paper.  The  uncertainty  coefficient 
represented  with  the  letter  "U"  is  a  normed  linear  measure 
of  the  analyst's  uncertainty.  If  there  is  no  uncertainty 
in  the  estimate,  U  =  0;  if  there  is  total  uncertainty  on 
the  whole  range  (i.e.,  all  values  are  equally  likely),  U  = 
1.  In .most  instances,  the  analyst  can  assign  a  reasonable 
value  to  U.  The  variance  of  x*  can  then  be  determined 
from  the  relationship: 

V (x*)  =  [(H  -  L)U]2/12,  0  <  U  <  1.  [6] 

It  is  difficult  for  the  analyst  to  visualize  the  dis¬ 
tribution  he  has  chosen  from  his  estimates  of  the  parame¬ 
ters.  However,  since  the  shape  of  the  beta  distribution 
is  determined  by  two  of  the  parameters  (a  and  b)  which  are 
in  turn  determined  by  the  values  of  the  mode  and  the  un¬ 
certainty  coefficient,  it  is  possible  to  get  a  reasonable 


1 0!  It  has  been  proposed  that  the.  analyst  assume  that  the 
r ange  of  the  cost  variable  1 S  equal  to  six  standard 
deviations ,  yielding  V { x* )  =  [(H  -  L)/6]2.  The  basis 
<5 o A.  this  assumption  Is  that  "most"  of  the  probability 
associated  with  distributions  such  as  the  normal  dis¬ 
tribution  Is  contained  In  the  Interval  ±  3  standard 
deviations  f rom  the  mean.  The  authors  of  this  paper 
do  not  Ind  this  to  be  a  very  strong  motivation.  The 
total  range  of  the  uniform  pdf  Is  contained  In  an  In¬ 
terval  of  ±  1.75  standard  deviations  and  this  distri¬ 
bution  Is  a  limiting  form  of  the  beta  distribution. 
Further,  the  significance  the  behavior  of  Infin¬ 
ite ,  symmetric  distributions  such  as  the  normal  pdf 
In  deriving  properties  of  the  [generally)  finite, 
asymmetric  beta  distribution  Is  questionable.  It 
seems  more  logical  to  allow  the  analyst  to  Input  his 
knowledge  of  the  uncertainty  via  the  uncertainty  co¬ 
efficient.  However,  If  the  user  prefers  this  device 
he  should  Input  the  value  U  «  0.577  when  using  pro¬ 
gram  SPET . 
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idea  of  the  distributional  form  from  a  set  of  normalized 
graphs  of  beta  pdfs  with  different  modes  and  uncertainty 
coefficients.  Such  a  set  appears  in  Appendix  B. 

Appendix  B  contains  ten  groups  of  graphs  of  normal¬ 
ized  beta  pdfs.  Each  group  contains  graphs  of  three  pdfs 
with  the  same  mode  but  distinct  uncertainty  coefficients. 
The  modal  values  vary  from  set  to  set,  beginning  with  .05 
and  increasing  by  .05  until  .50  when  reading  the  abscissa 
from  left  to  right,  or  beginning  with  .95  and  decreasing 
by  .05  until  .50  when  reading  from  right  to  left.  To  use 
Appendix  B,  the  analyst  computes  the  estimated  normalized 
mode  (N)  from  his  estimates  of  L,  H,  and  M  with  the  rela¬ 
tionship 


N  "=  (M  -  L)/(H  -  L) 


[7] 


and  then  selects  the  group  of  pdfs  whose  normalized  mode 
is  closest  to  this  computed  N.  From  the  three  graphs  in 
this  subset  the  analyst  can  see  how  the  pdf  varies  with 
the  uncertainty  coefficient  and  get  a  reasonable  idea  of 
the  shape  of  the  distribution  he  has  chosen.  Alternative¬ 
ly/  he  may  look  at  the  set  before  choosing  the  uncertainty 
coefficient  and  use  the  information  he  gains  to  help  him 
select  the  value  for  U. 

Example 

Assume  an  analyst  is  studying  the  cost  of  some  sub¬ 
system  "S".  By  analogy  with  other  systems,  or  by  some 
other  technique,  he  determines  that  the  cost  of  S  will  be 
something  greater  than  $7,000  but  less  than  $12,000.  Fur¬ 
ther,  utilizing  a  CER,  or  by  some  other  technique,  he  es¬ 
timates  its  most  likely  value  at  $10,500.  He  calculates 
the  normalized  mode  "N"  from  the  equation: 


M  -  L  10,500  -  7,000 
=  H  -  L  =  12,500  -  7,000  ~ 


The  set  of  graphs  corresponding  to  this  system  is  found 
between  pages  B-19  and  B-21.  These  graphs  are  read  from 
right  to  left. 

After  the  analyst  has  developed  the  distributions  of 
each  subsystem  he  faces  the  problem  of  determining  the 
distribution  of  their  sum  (the  total  cost  of  the  system) . 


11 


-  -  -  » 


yrr 

-  - *  ; 


*sr 


■  ...  . 


o 


o 


COMBINING  SUBJECTIVE  PROBABILITY  DISTRIBUTIONS 

Of  the  several  techniques  that  have  been  employed  by 
cost  analysts  to  combine  statistical  distributions  repre¬ 
senting  their  subjective  probability  distributions,  two  of 
the  more  popular  are  derivation  of  moments 11/  and  Monte 
Carlo  simulation. 12/  Each  of  these  techniques  has  its  ad¬ 
vantages:  the  former  can  be  done  with  tables  and  a  desk 

calculator,  whereas  the  latter  requires  access  to  an  elec¬ 
tronic  computer.  The  latter,  however,  is  much  faster  and 
easier  to  use.  For  this  reason  Monte  Carlo  simulation  is 
the  technique  employed  in  this  research.  Appendix  C  docu¬ 
ments  Program  SPET,  a  computer  program  for  adding  indepen¬ 
dent  statistical  distributions  by  means  of  Monte  Carlo 
simulation  (SPET  also  performs  other  calculations  dis¬ 
cussed  in  the  next  section) . 

Program  SPET  has  been  used  successfully  on  an  inter¬ 
active  time-sharing  computer  system.  Basically  the  user 
enters  the  parameters  of  the  statistical  distributions  se¬ 
lected  by  the  analyst  and  the  program  generates  frequency 
distributions  and  summary  statistics  of  the  total  system 
cost.  Details  and  an  example  of  the  inputs,  outputs,  and 
operation  of  the  program  can  be  found  in  Appendix  C. 

THE  INDEPENDENCE  ASSUMPTION:  A  PROBLEM 

When  the  analyst  decomposes  the  system  he  is  costing 
into  several  subsystems,  it  is  very  unlikely  that  the 
costs  of  the  various  subsystems  are  always  statistically 
independent  of  one  another.  For  example,  the  cost  of  the 
propulsion  system  of  a  guided  missile  is  probably  corre¬ 
lated  with  the  cost  of  its  payload.  The  correlation  would 
be  positive  if  an  increase  in  payload  cost  was  due  to  an 
increase  in  payload  size  which,  in  turn,  would  require  a 
more  powerful  and  hence  more  costly  propulsion  system.  On 
the  other  hand,  the  correlation  would  be  negative  if  an 
increase  in  payload  cost  was  due  to  a  reduction  in  payload 


1 1  /  See  S.  Sobel,  A  ComputeAized  Technique,  to  ExpAe&i  Un- 
ccAtainty  in  Advanced  System  Co&t  E itimate*  [Bed^oAd, 
JffT:  The TO  CoJLp.,  TW-372  & ,  196  3),  and  W.  H.  Suth- 
e Aland,  A  Method  {oa  Combining  A &ymmetAic  ThAee-Value 
PAedictionA  oj  Time  o a  Co6t  (McLean,  iTXl  R e&eaAch 
A nalyiii  CoAp.,  KAC-P-65,  T 972). 

12/  P.  F.  Pienemann,  ojo.  c it,  and  P.  F.  Schoe^eA,  op. 
c  it. 
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size  brought  about  by  miniaturization  which,  in  turn, 
would  require  a  less  powerful  and  hence  less  costly  pro¬ 
pulsion  system.  It  is  impossible  to  determine  a  priori 
whether  the  correlation  between  the  costs  of  any  two  sub¬ 
systems  is  positive  or  negative,  but  the  experience  of  the 
authors  suggests  that  in  the  majority  of  cases  it  will  be 
positive. 

If  an  analyst  assumes  that  the  statistical  distribu¬ 
tions  (random  variables)  representing  subsystem  costs  are 
independent  when  in  fact  they  are  positively  correlated, 
he  will  underestimate  the  variance  of  their  sum.  To  see 
this,  consider  the  expression  for  the  sum  (S)  of  the  vari¬ 
ance  of  n  random  variables  (X^) : 

n  n-1  n 

V(S)  =  ZVUJ  +  2  E  Z  CovfXi,  X.)  18] 

i=l  i=l  j=i+l  3 


The  assumption  of  independence  implies  that  Cov(X^,  X^)  = 
0  for  all  which  in  turn  implies  that  the  expression 


n-1  n 

2  E  E  Cov(Xi#X.»)  =  0. 
1=1  j=i+l  J 


IF  the  Xi  are  dependent,  this  term  could  be  positive,  neg¬ 
ative,  or  zero.  If  it  is  positive,  its  deletion  from  [8] 
by ^assuming  independence  results  in  an  understatement  of 
V(S).  An  understatement  of  V(S)  could  be  a  serious  prob¬ 
lem  because  it  results  in  a  confidence  interval  about  the 
mean  of  the  total  system  cost  distribution  that  is  smaller 
than  it  should  be.  This  might  cause  the  decision  maker  to 
posit  unwarranted  confidence  in  the  estimate. 


Assuming  the  random  variables  representing  subsystem 
costs  are  positively  correlated,  the  magnitude  of  the  un¬ 
derestimate  of  V(S)  is  directly  related  to  the  number  of 
variables  in  the  sum.  To  illustrate  this  fact,  consider 
the  following  example: 


There  are  two  positively  correlated  random  variables, 
X^  and  X2.  Then 

V(S)  =  V(XX)  +  V(X2)  +  2Cov(X1,X2)  19] 

Now  assume  X^  *  +  Z2. 


Then 


O 


c 


V(S)  =  V(Zj)  +  v(z2)  +  V(X2)  +  2  CovIZ^Zj)  + 

-  -  2  Cov(Zi,X2)  +  2  Cov(Z2,X2)  -  - - 

-  V(ZX)  +  V{z2)  +  V(X2)  +  2  Cov(X1,X2)  + 

'  -  -  -  2  Cov(Z1,Z2)  110] 

Assuming  independence  in  the  case  of  two  variables 
results  in  the  deletion  of  2  Cov(Xi,X2)  from  V(S) . 

-  -  However,  in  the  case  of  three  variables  the  assump¬ 
tion  of  independence  results  in  the  deletion  of  not 
;  ;  only  2  Cov(X^,X2)  but  2  Cov(Zi,X2)  as  well.  Clearly, 

■  if  the  covariances  are  positive,  V(S)  is  understated 
more  in  the  case  of  three  variables  than  in  the  case 
of  two. 

Therefore  in  the  case  where  the  random  variables  repre¬ 
senting  subsystem  costs  are  positively  correlated  not  only 
is  the  variance  of  the  total  system  cost  underestimated, 
but  the  magnitude  of  the  underestimate  is  directly  related 
to  the  number  of  variables  making  up  the  sum. 

The  obvious  solution  to  this  problem  is  to  incorpo¬ 
rate  into  computer  Program  SPET  provisions  for  the  consid¬ 
eration  of  probable  correlation  among  the  variables  and 
then  proceed  to  estimate  the  nature  of  the  correlation. 
Although  the  former  idea  presents  no  problem,  the  latter 
appears  to  be  a  most  difficult  task.  It  is  not  clear  at 
this  point  how  to  systematically  estimate  the  correlation 
among  the  variables  representing  subsystem  costs.  Hope¬ 
fully  a  credible  technique  for  doing  such  will  become  ap¬ 
parent  to  someone. 

Program  SPET  has  been  designed  to  provide  some  in¬ 
sight  into  the  significance  of  the  independence  assumption 
in  two  ways.  First,  by  using  the  same  random  number  in 
sampling  from  all  the  subsystem  distributions,  a  distribu¬ 
tion  of  total  cost  which  the  printout  titles  "Dependent 
Beta"  is  derived.  The  technique  imposes  a  functional  re¬ 
lationship  between  all  the  variables.  Note  that  this 
functional  relationship  is  not  an  arbitrary  relationship 
but  is  imposed  by  the  forms  of  the  subsystem  pdfs  devel¬ 
oped  by  the  analyst.  Specifically,  if  Fi(Xi)  is  the  cumu¬ 
lative  distribution  function  of  the  ith  subsystem  then: 

n  n 

E  X.  =  X.  +  E  Fx-MF.fX.H  [11] 

i=l  1  1  i=2  1  1  1 
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where  Fi-*-  is  the  functional  inverse  of  Ff.  The  functions 
Fi  and  Fi-1  are  much  too  complex  to  derive  the  specific 
form  of  the  relationship  imposed  but  the  technique  does 
impose  a  very  real  positive  correlation  between  the  varia¬ 
bles. 


As  a  second  means  of  examining  uncertainty  without 
imposing  the  independence  assumption,  SPET  prints  the  sta¬ 
tistics  for  a  uniform  distribution  of  total  cost  on  the 
interval  between  the  sum  of  the  minimum  costs  of  the  sub¬ 
systems  and  the  sum  of  the  maximum  cost  of  the  subsystems. 
This  is  meant  to  serve  as  a  "worst  case."  SPET  also 
prints  the  total  cost  distribution  assuming  each  of  the 
subsystem's  costs  is  uniformly  and  independently  distrib¬ 
uted. 


SUMMARY 


them: 


Several  conceptual  points  have  been  discussed,  among 


(1)  the  nature  of  uncertainty  in  cost  analysis; 

(2)  the  value  of  treating  uncertainty  explicity  in  cost 
analysis; 

(3)  a  problem  inherent  in  using  the  classical  linear  re¬ 
gression  model  as  a  basis  for  statements  on  cost  un¬ 
certainty; 

(4)  the  subjective  nature  of  cost  uncertainty; 

(5)  t'ne  properties  of  the  beta  distribution  and  how  they 
can  be  used  to  facilitate  cost  uncertainty  analysis; 
and 

(6)  the  dependence  of  subsystem  costs  and  its  impact  on 
statements  about  uncertainty. 

• 

The  basic  practical  contribution  of  this  paper  is  a 
computer  program  for  generating  statements  on  cost  uncer¬ 
tainty  that  permits  the  analyst  to  input  any  imaginable 
probability  distribution  to  represent  a  subsystem  cost. 


*  fl 


APPENDIX  A 


ESTIMATION  OF  PARAMETERS  a  AND  b 
COMPUTER  PROGRAM  PARAM 


o 


Introduction 


One  way  an  analyst  can  estimate  parameters  a  and  b  of 
the  generalized  beta  probability  density  function  (pdf)!/ 
is  to  simultaneously  solve  two  functionally  independent 
equations  in  a  and  b.  Consider  the  mode  —  the  "most 
likely"  value  —  of  the  generalized  beta  pdf: 

M(x*)  =  (aH  +  bL)  /  (a  +  b)2/  [Al] 

» 

Using  a  cost-estimating  relationship,  or  other  methods, 
the  analyst  can  estimate  the  mode  (M)  of  the  subsystem 
cost  he  wants  to  represent  with  a  beta  pdf.  Substitution 
of  M  into  equation  [Al] ,  along  with  L  and  H,  yields  one 
equation  in  a  and  b : 

M  '=  (aH  +  bL)  /  (a  +  b)  [A7] 

Another  equation  in  a  and  b  can  be  obtained  from  an 
estimate  of  the  probability  (P)  that  the  cost  of  the  sub- 


y  See  equation  [2]  on  page.  8. 

y  Thlt  exprettlon  can  be  obtained  by  solving  dlglx.*)]/ 
dx*  *  0  v*,  where  git*)  It  given  by  equ alien  [2]. 

However,  It  It  eatler  to  obtain  exprettlont  £ ok  the 
mode  at  u)zll  at  the  mean,  E(x*),  and  the.  variance, 
Vlx*),  0(5  the  generalized  beta  pdfi  by  meant  o{  tlmple 
algebraic  operatlont  on  the  exprettlont  &or  thete  tta- 
tlttlct  derived  firom  the  mo> te  £ amlllar  normalized  beta 
pd£.  The  mode,  mean,  and  variance  o<J  the  normalized 
beta  pd£  are: 

M ( x )  =  a/[a  +  b)  [A2] 

E ( x)  «  (a  +  J)/la  *  b  +  2)  (A3] 

V[r)  «  [la  +  i)(b  +  I )  3  /  f  C  a  +  b  *  2 ) 2  (a  +  b  +  3)]  (A4) 

Note  that  M(x)  *  a/ (a  +  b)  *  M[(x*  -  L)/(H  -  L ) ]  - 
I M ( x * )  -  l]/ [H  -  L) ;  tolvlng  ior  M(x*)  yleldt  equation 
[Al] .  Proceeding  tlmllarly  $or  E  ( x* )  and  1/  ( x* )  yleldt 

E(k*)  •  [aH  +  bL  *  H  +  l)/(a  +  b  *  2)  [A5J 

V(x*)  -  [(a  +  mb  *  1)IH  -  L) z]  /  I(a  +  b  ♦  2 ) 2  (a  +  b  + 
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system  will  lie  within  a  subinterval  of  its  range  [L,  H]  . 
For  convenience,  this  subinterval  is  taken  as  the  interval 
from  L  to  the  midpoint  between  L  and  M,  i.e.,  [  (L  +  M)/2] . 
Then  an  equation  in  a  and  b  results  from  the  relationship 

L+M 

P  =  IC/(H  -  L)a+b+1]  r*~  (x*  -  L)a  (H  -  x* ) b  dx*  [A8] 

L 

# 

where  C,  x,  a,  and  b  are  defined  in  equation  [2] . 

Equations  [A7]  and  [A8]  comprise  two  functionally  in¬ 
dependent  equations  in  a  and  b,  and  when  solved  simultane¬ 
ously  determine  unique  values  (estimates)  for  a  and  b. 


Program  Description 

Program  PARAM  is  written  in  FORTRAN  IV  for  use  in 
conjunction  with  a  PDP-10  computer  in  an  interactive  time¬ 
sharing  mode.  It  is  designed  to  solve  equations  [A7]  and 
[A8]  for  parameters  a  and  b,  given  values  for  L,  H,  M,  and 
P.  This  is  accomplished  in  the  following  sequence: 

1  -  The  inputs  L,  H,  and  M  are  normalized.  Denoting 
the  normalized  counterparts  of  L,  H,  an-  M  by  1-  h.-  and  m, 
respectively, 

1  =  0 

h  =  1 

m  =  (M  -  L)/(H  -  L) 

2  -  A  standard  root-finding  technique—/  is  employed 
in  Subroutine  Beta  to  find  the  \alues  for  a  and  b  that 
satisfy 


m/2 

P  -  C  ./  xa(l  -  x)b  dx  [A9] 

n 

where  P  is  supplied  by  the  user,  m  is  determined  from  a 
user  supplied  datum  (M) ,  and  c,  x,  a,  and  b  are  as  defined 
in  equation  [2] .  Note  that  [A9J  is  the  normalized  version 
of  [A8] . 


3/  The  Koot-  finding  -technique  iA  known  <u  the.  n^aJtAZ  po- 
Aition”  method  and  La  iound  in  moAt  textA  on  nu.me.Ki cal 
analyAiA. 


A-2 
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Subroutine  Beta  calls  Function  Subprogram  Gamma  to 
compute  values  for  T(n).  This  is  accomplished  by  use  of 
the  relation 

T(n+1)  =  nr (n)  [AlO] 

and  an  interpolation  procedure. 

By  repeated  application  of  [AlO],  r  (n)  for  any  n>0 
can  be  expressed  as  the  product 

(n  -  1) (n  -  2) . ..T(n*)  [All] 

where  1  -  n*  -  2.  Since  T(n)  is  well-behaved  in  the  range 
1  5  n*  -  2,  the  values  of  T(n*)  can  be  approximated  using 
a  "table  look-up"  interpolation  procedure.  Function  Sub¬ 
program  Gamma  uses  such  an  interpolation  device  in  con¬ 
junction  with  relation  [All]  to  evaluate  T  (n) . 

3  -  The  four-tuple  (a,b,L,H)  is  printed  as  output. 


User  Instructions 


The  following  example  demonstrates  the  use  of  Program 
PARAM.  Note  that  the  user  supplied  portions  of  the  exam¬ 
ple  are  underlined: 


•  «Ji\  PArtAo 


-  cSflMAi'cS  Or  L,H,MtP 
JtJfc;,  I4fc :0,  lOoO'O..  14 


ALPHA  IS  1.07 

i3  Lf  A  IS  I .  D0 

L0«  Id  dM-J.CQ 
ill  On  IS  14010.00 


OP  UiMl  I  d  2 

f 

lXH 


A- 3 


e 


Program  Listing 


B£dl2 

Md2*) 

J0d  33 

30d4J 

30d  53 

00d  60 

3387  0 

30380 

00a  90 

00900 

009  1 0 

00920 

00930 

00940 

00950 

00y60 

00970 

009d0 

00993 

31000 

«  I  *  %  |  .  a 

•  O  I  V 

31320 
31330 
31340 
31053 
31360 
31370 
3  1283 
31390 
31100 
01.1  13 
31  120 
01130 
3.1143 
0.1 1  50 
011 60 
31173 
01180 
31190 


A=BU 

220  PP=0. 

MM  1 1  U=0. 

HLl'JKN 

50  A’RIJtUI  .60) 

63  rOfiViAl  (  Ih-  ,28H3ETA  PARAMETERS  OUT  OF  RANGE) 

RETURN 

end 

c 

FUNCTION  GAMMA (GP) 

REAL  G ( 22 ) 

DATA  G/l . ,.97 35,. 95 135,. 93304,. 91 8 1  7, 

-+  .  9064  , .  897 47  ,  . 8  9  1 1 5  , .  88725  , .  38565 , 

-+  . 88623  , .  88 887  ,  .893  5  2  ,.90012,.  90864 , 

+  .51  936,.  93138,.  94  56  I  ,.96  177,.  97  933,  1  .  ,  1  ./ 

11=16 
JJ  =  16 
ERROR =3. 

SU1A=t. 

IF (GP-57 .4 )  10,10,20 
13  IF vGr- 1 .F0C-20)  23,30,30 

20  MriIT£(  1 1  ',43 ) 

40  FORMAT (  IH-  ,2d.HGAMMA  PARAMETER  OUT  OF  RANGE) 

ERH0R=cRk0R4  I  . 

RETURN 

30  IF (GP-2. )  53,50,60 

63  SUM  =  SUM* ( GP - 1 . ) 

GP=GP- 1  . 

GO  TO  30 

50  THGP-I.)  70,33,80 
70  SUM=SUM/GP 

GP=GP(- 1 . 

80  1  =  ( GP-  I .  )/  .  05  +  1  . 

XI =1-1 

CPL=! .+XI *.  05 

CAMFN=G(1 )+(G( I+l )-C(IJ )*<GP-CPL)/ .05 

GAMMA  =OA.V.F,\'*SUM 

RETURN 

END 


A-4 


304  1  0 

IF (GP-o7 . 4 )  43,40,50 

03420 

40 

Cl  =GA.VMA  (GP) 

304  30 

UP=A+ 1 . 

00440 

C2=GAf.;.'A  (GP ) 

004  50 

C3  =  CI/C2 

00460 

CP=d  + 1 . 

004  70 

C4=3A,V.V.A(GP) 

00460 

C=C3/C4 

00490 

70 

X/.,M=XM.Vi+2. 

00500 

D)  80  1=2,3 

005  1-0 

FT(  I  )=f  (I.)**A*(  l.-TCI  ) )  **3 

0G520 

T1I+1  )=T(  I  >  -jSP 

00530 

<30 

CONTINUE 

00540 

1(2)  =1(4) 

00550 

XSUM=FT ( 1 )+4.*FT(2)+FT(3) 

00560 

Fid  )=rT ( 3) 

00570 

XI  NT  LG=XI NTEG+XSUM 

003  00 

IF  (  XMM-I  NT  EG)  70,90,90 

005  90 

90 

AA=P-C*5P/3.*XINTEG 

006  £0 

IF ( AdS( AA) - . 0?3 1 )  100,100,1 

026  !  0 

110 

IF ( AA-AAS)  120,130,133 

000  20 

130 

IF(AA)  140,143,150 

006  30 

140 

AA  1  =AA 

00640 

GP  1  =B 

00o50 

GO  TO  160 

00660 

150 

AA2=AA 

006  7  0 

GP2  =B 

006  8  0 

160 

GO  TO  (  170, 180)  ,JJI 

006  90 

120 

JJ  1  =2 

007  00 

GO  TO  130 

007  1  0 

170 

B=W 

00720 

rt=w* 2. 

007  30 

GO  TO  190 

007  40 

160 

3-( AA2*GF 1- AA1 *CP2 )/ ( AA2-AA 1 

007  50 

193 

AASMA 

00760 

GO  TO  30 

007  70 

1  03 

IF ( PP-NMl DD)  200,200,210 

00780 

200 

CO  TO  220 

007  93 

210 

3B=B 

006  00 

B-A 

o 


003  1  0 

RbAL  L,M 

30J23 

11  =  16 

30 

JJ  =  16 

02040 

writ  bur,  i0> 

020  5  A 

10 

FOHMAK  1H-, 20.-1  EbT IV. AThS  OF  L,H,M,P/> 

02360 

RLAD(  JJ  ,20 )  L,H,M,P 

00070 

20 

FORMAT  C4F) 

02060 

CALL  bLTACL.H.V  ,P,A,B> 

003  V0 

WftITbd  I  ,30)  A,B  ,L,H 

0a  1  00 

30 

FORMAT ( 1H-,9HALPHA  IS  .F8.2/I0H  BETA  IS  ,F8.2/ 

00110 

+  1BH  LOrt  IS  .F8.2/10H  HIGH  IS  ,F8.2) 

00120 

rtHiTEI I I ,407 

00130 

40 

FORMAT ( 1H-) 

00140 

LiND 

00150 

C 

00160 

SUBROUTINE  BETA (L,H , M,P ,A ,B) 

00  1  70 

R  LAL  L  ,M ,  I  NT  LG  ,  NV.ODb ,  N,M I  D ,  NMI DD  ,  FT  ( 3 )  ,  T  ( 4 ) 

00  1  80 

11=16 

22 1  V0 

JJ  =  16  J 

002  00 

INI  LG =4G. 

002  1 3 

S  LARCH  =22 . 

00223 

NMOCb=(.v,-L)/  (H-L) 

00233 

NMI  D=NM0Dc/2. 

00240 

IMP-NMID  10, 10,20 

002  50 

20 

PP=P 

00260 

NMI rD=NMl D 

00270 

•P=l.-P 

00230 

NMI  0=1  .-NMI D 

002  90 

NV.()DE=  I  .-NMOrE 

003  20 

10 

SP=NMI D/I  NT EG 

00310 

AAS=-l . 

00o20 

JJ  1  =  1 

00330 

'A’=S  LAkCH*  ( 1  ,-NV.ODE) 

00340 

6=0. 

00353 

30 

A=B*N.V,0DT/  (  1  .-NMOOE) 

0eiO0 

XMM=3. 

00370 

XI  NTcG=0. 

023  d0 

1(2 ) =SP 

003  VO 

f  T  ( 1 ) =3. 

004  23 

GP=A-kJ  +2 . 

I 
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APPENDIX  B 


GRAPHICAL  AIDS  FOR  SELECTING 
THE  UNCERTAINTY  COEFFICIENT 


-TO- RIGHT  HIGH  VARIANCE  RIGHT-TO-LEFT 

alizecl  Mode  =  .05  Uncertainty  Coefficient  =  .75  Normalized  Mode 


TO-RIGHT  LOW  VARIANQ5  RIGHT-TO-LEFT 

lized  Mode  *  .10  Uncertainty  "Coefficient  -  .25  Normalized  Mode 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


TO- RIGHT  MEDIUM  VARIANCE  RIGHT-TO-LE 

lized  Mode  »  .15  Uncertainty  Coefficient  »  .50  Normalized 


xxxxxxxxxxxx 


•  LEFT-TO- RIGHT  MEDIUM  VARIANCE  RIGHT-TO-LEFT 

Normalized  Mode  «*  .20  '  Uncertainty  Coefficient  ■  .50  Normalized  Mo 


xxyxxxxxxxxxxxxyxxxxxxxxx> 


-fcgjgg  j#i 


xxxxxxxxxxxxxxxxxxxxxxx 


TO- RIGHT  HIGH  VARIANCE  RIGHT-TO-LEFT 

lized  Mode  =»  .45  Uncertainty  Coefficient  =  .75  Normalized  Mo* 


in 

in 
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<u 
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X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


LEFT-TO- RIGHT  MEDIUM  VARIANCE  RIGHT-TO-LEFT 

Normalized  Mode  =  .50  Uncertainty  Coefficient  *  .50  Normalized  Mode 
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Descrip  >un 


Program  SPET  is  written  in  FORTRAN  IV  for  use  in  con¬ 
junction  with  a  Univac  1108  computer  in  an  interactive 
time-sharing  mode.  It  is  designed  to  approximate  the  fre¬ 
quency  distribution  of  the  sum  of  up  to  fifty  independent 
beta  or  other  random  variables  by  means  of  Monte  Carlo 
simulation.  This  is  accomplished  in  the  following  manner: 

» 

(1)  The  user  specifies  the  beta  variables  by  entering 
the  lowest  (L) ,  most  likely  (M) ,  and  highest  (H)  value  as 
well  as  the  uncertainty  coefficient  (U)  of  each  variable. 
He  specifies  other  variables  by  entering  the  lowest,  most 
likely,  and  highest  values,  as  well  as  the  mean,  variance, 
and  cumulative  density  function  (cdf)  of  each  variable. 

He  enters  only  those  values  of  the  cdf  F(x)  corresponding 
to  x  =  L  +  i (H  -  L)/10,  i  =  1,...,10. 

(2)  For  all  beta  variables  the  computer: 

(a)  Computes  beta  parameters  a  and  b  from  M  and  U 
by  first  converting  U  to  V  with  equation  [6]  and  then  si¬ 
multaneously  solving  equations  [3]  and  [4] . 

(b)  Computes  the  discrete  cdf,  F(x),  for  x  = 
i/10,  i  =  1,...,10,  using  Subroutine  DQG32 . W 

(3)  Next  the  class  intervals  for  the  distribution  of 
the  sum  of  the  beta  and  other  variables  are  computed. 
Adding  the  L  values  and  throughput  (TPUT)l/  establishes 
the  lower  limit  of  the  range;  adding  the  H  values  and  TPUT 
establishes  the  upper  limit.  The  range  is  then  divided 
into  15  intervals  of  equal  width.—' 

(4)  Four  frequency  distributions  of  the  sum  are  gen¬ 
erated.  The  first  distribution  results  from  the  assump¬ 
tion  that  the  distributions  making  up  the  sum  are  statis- 


1J  VQG-32  ui>  £6  the  32  point  GauAAian  qua.dna.tu.ne.  method  o£ 
integration.  It  iA  taken  {Horn  Convolution  o<(  JnverAe 
Beta  ViAtributionA  by  a  Sampling  Technique  [BetheAda, 
Ml) :  Mathematica,  1  n  c.  TTYT ). 

2J  Throughput  meanA  conAtant,  and  iA  usually  the  coAt  o& 
a  AubAyAtem  that  iA  known  with  certainty. 

3J  The  number  o {  intervalA  can  be  varied  by  aAAigning  the 
deAired  number  to  KK  in  line  11  oi  the  main  program. 
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tically  independent.  It  is  generated  as  follows: 

(a)  Obtain  a  random  number  lying  between  zero  and 
one  from  a  uniform  random  number  generator.!/ 

(b)  Compare  the  generated  number  with  the  values 

of  the  discrete  cdf,  F(x),  for  one  of  the  variables  and 
note  the  interval  [F  (x^)  ,F  (x-j )  ] ,  <  xa  ,  into  which  it 

falls. 


(c)  Find  the  x-value,  xj_  <  x  <  x j ,  corresponding 
to  the  random  number  by  means  of  linear  interpolation.^,/ 

(d)  Transform  the  x-value  from  its  normalized 
(0,1)  value  to  its  standard  (L,H)  value  by  means  of  the 
transformation  x*  =  L  +  (H  -  L)x.  [Cl] 


4_/  An  uniucc.e.64)  hut  attempt  wa6  made  to  fiind  a  machine- 
independent  nandom  numben  genenaton  Ion  inclusion  in 
the  pnognam.  Thene&one,  Pnognam  SPET  tiegu.in.ei>  the  use 
0(5  a  ui  en- supplied  genenaton.  Make  the  appnopniate 
change-i  in  line  150  0( 5  the  main  pnognam  to  accommodate 
the  genenaton  (it  may  be  neceA&any  to  al&o  change 
Unci,  68-69,  155,  158,  162,  175,  1  78,  18 2,  262,  and 
370-372 ] . 

5/  Siep-S  (fa)  and  (c)  can  be  illuttnated  fion  a  nonmalized 
nandom  vaniable  x  a-6  ^ollom: 


F  (x) 
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(e)  Repeat  steps  (a)  -  (d)  for  every  variable. 

(f)  Compute  the  observation, 

N 

X.,  =  TPUT  +  Z  Xt 
J  i  —1  A 


where  N  denotes  the  number  of  random  varia¬ 
bles. 

(g)  Find  the  class  interval  [see  (3)  above]  in 
which  X-;  occurs  and  register  one  occurrence  in  that  inter¬ 
val. 


(h)  Repeat  steps  (a)  -  (g)  for  j  =  1,...,M  obser¬ 
vations,  where  M  =  number  of  desired  observations.  When 
step  (h)  is  terminated,  one  has  a  frequency  distribution 
of  the  sum  of  independent  random  variables. 

(5)  The  second  distribution  of  the  sum  is  generated 
in  the  same  way  as  the  first  except  for  step  (e) ,  which 
should  now  read,  "Repeat  steps  (b)  -  (d)  for  every  varia¬ 
ble.  "  This  means  that  the  same  random  number  is  used  to 
obtain  an  observation  on  each  component  of  the  sum  rather 
than  a  new  number  as  done  when  constructing  the  first  dis¬ 
tribution,  The  procedure  of  using  only  cnc  random  number 
introduces  a  correlation  among  the  component  variables  be¬ 
cause  when  the  value  of  one  variable  is  known,  it  in  turn 
maps  uniquely  to  the  values  of  all  other  variables.  This 
correlation  is  positive  because  the  cdfs  are  positive 
mono tonic  functions. 

(6)  The  third  distribution  is  simply  the  uniform  dis¬ 
tribution  over  the  range  of  the  first  distribution.  Both 
the  second  and  third  distributions  serve  as  indicators  of 
how  a  violation  of  the  independence  assumption  could  af¬ 
fect  the  distribution  of  the  sum  and  the  summary  statis¬ 
tics. 


(7)  The  fourth  distribution  is  the  same  as  the  first 
distribution  except  the  component  variables  are  all  uni¬ 
form  random  variables.  This  distribution  serves  as  an  in¬ 
dicator  of  the  relative  sensitivity  of  the  distribution  of 
the  sum  to  the  degree  of  uncertainty  in  the  component  var¬ 
iables. 

(8)  Along  with  the  four  frequency  distributions  just 
described,  the  computer  generates  the  mean,  mode,  standard 
deviation,  variance,  90%  confidence  interval  about  the 
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mean,  and  the  probability  of  exceeding  the  mean  for  each 
of  the  distributions.  In  addition,  the  user  can  specify 
any  confidence  interval  about  the  mean  and  the  probability 
of  exceeding  any  number  within  the  range  of  the  distribu¬ 
tion  and  the  computer  will  generate  it  for  him. 

User  Instructions 


The  following  example  illustrates  the  use  of  Program  SPET: 


«SPtf 


mix 

••example** 

tow* lr  or  uisruiuurioNS 

* 

DATA  FItem  OR  TtKMlNALiai  INPUT 

2 

LOWEST » MOIH.'  *  HIGHEST  •UNCERTAINTY  COfcFF 

DATA 

100*2SO*'300»  .8 
DATA 

500o3r»0»590o.2 

OATA 

200»400»67b»  ,5 
OATA 

Ib0»3>0*tf00# .4 

NUMBER  Of  OTmER  DISTRIBUTIONS 

1 

LOaEST » M00C  o HIGHEST o WEAN* VARIANCE* DISCRETE  COF 

DATA 

50# 75. 100.  7b.  208.5#  .!•  *2#.3»  ,4.  .5*  .6.  .7.  .0.  .9. 1.0 

THROUG'IPUT 

55 

NU'-'iJEK  Of  0  JSERVA.  IONS 
20D0 

SEED  RANDOM  NUMBER  GENERATOR 
87654521 


••EXAMPLE »» 


INPUTS 

OBSERVATIONS  =  2000 

SEED  =  87654521 


BETA  VOLE 

LOWEST 

mooc 

HIGHEST 

U  CCEFf 

NMODE 

alpha 

BETA 

MEAN 

VARIANCE 

1 

100.0 

250.0 

900.0 

•  60 

.19 

.3 

l.J 

396.* 

33006.0 

2 

500.0 

350.0 

390.0 

•  20 

.56 

39.6 

31.6 

3*9.9 

27.0 

3 

200.0 

400.0 

675.0 

•  50 

.42 

3.7 

5.1 

*06.9 

*697.* 

A 

150.0 

330.0 

eoo.o 

.**0 

.20 

3.6 

9.4 

3*9.3 

561*. 6 

OT»CR  VHLE 

s 

50.0 

75.0 

100.0 

75.0 

206.3 

THROUGHPUT 

53.0 

35.0 

35.0 

••SUMS** 

055.0 

2900.0 

1612.5 

**353.2 

o?»cn  vBtx  ‘  cor 

%  *1000  .2000  .5000  .4000  .5000  .6000  .7000  .0000  .9000  1.0000 
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INTERVAL 


OUTPUT 


INDEPCNOCNT  OETA/OTHER  DEPENDENT  IIETA/OTHER  TOTAL  UNICORN 


835.0  - 
972. »  - 
1J)0.3  - 
12X8.0  - 

1303.7  - 

1573.3  - 

1661.0  - 

1708.7  - 

1930.3  - 
2071.0  - 

2711.7  - 
23x9.3 
2X87.0  - 
2A2X.7  - 

2762.3  - 


972.7 

1110.3 
12X8.0 

1385.7 

1523.3 
1661 .0 

1798.7 

1936.3 
207X.0 

271 1.7 
23X9.3 
7X87,0 
267X.7 

2762.3 
2900.0 


TOC  COC 


0. 

.onto 

.02X0 

.11X0 

.2770 

.2300 

.19X5 

.12X0 

.060'. 

.0185 

.0010 

0. 

0.  ’ 

0. 

0. 


.0010 

.0250 

.1390 

.3600 

.6020 

.7965 

.9205 

.9805 

m  QxiuO 

l.or.  jc 
1 .0000 
1.0013 
1.0000 
1.0000 


*  .0110 
.05X5 
.10X0 
.12X5 
.1350' 
.1365 
.1170 
.11)5 
.0875 
.05X0 
.0385 
.0175 
.0065 
.0015 
0. 


COF 

POF 

COF 

POF 

•  0110 

.0667 

.0667 

.0010 

0  0655 

.0067 

.1333 

.0025 

.0667 

.2000 

•  0215 

,?940 

.0667 

.2667 

.0X35 

•  429b 

.0667 

.3333 

,0805 

.0660 

.0667 

.  X  00G 

.11X0 

.6630 

,0667 

.46**7 

.  1575 

.794  5 

.06u7 

.5333 

•  16]  5 

.0667 

.60U0 

.  1525 

.9360 

.0667 

.6**67 

.1070 

.9745 

.0667 

.7333 

.0395 

•9‘J20 

.0667 

.8000 

.0X15 

.9**115 

•  06*67 

.8667 

.0210 

•  nouo 

.06**7 

.9333 

.03 60 

•  0000 

.0667 

1.0000 

.0305 

INDEPENDENT  IAUFOx 


.0010 
.0035 
.0250 
.0685 
.1X90 
.2630 
,X2J>5 
.5820 
.73X5 
.8x15 
.9319 
.9725 
.9:35 
.9  >"S 

1.030: 


..KEAN** 

..MODE** 

..VARIANCE*. 

•«STO  DEVIATION*. 

*.901  CONFIDENCE  lNlERVAL*- 


1615.5 

1592.2 

4X8X3, 7 

•  211.8 

1307.2. 
,  o 


1621.1 

1592.2 

128X88.0 

358.5 

J12X.B. 

236X.1 


355352.1 

596.1 


938.3. 

3796,8 


1869.2 

1867.5 

103015.0 

321.0 

1337.9, 

?«23.» 


..PROS  EXCEED  MEAN.. 


ANOTHER  CONFIDENCE  INTERVAL? 
«?95*  CONFIDENCE  INTERVAL** 


1277.0. 

2218.5 


1071.2. 
2900.0  . 


886.6. 

28X9. X 


1258.3. 

252X.7 


another  confidence  interval? 

0 


prob  c retro  some  value? 

1450 

..PROU  EXCEEO  1X50.0** 


PROD  EXCEED  SOME  VAL1*? 

0 


additional  observations? 

0 


i 


anothcr  seed? 

• 

STOP 

SRU'SIIX.7 


.  *4-  *». 
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(1)  The  user  enters  a  title  up  to  60  characters  long. 

(2)  If  he  is  using  beta  distributions  to  represent 
some  or  all  of  his  component  distributions,  the  user  en¬ 
ters  the  number  of  distributions  to  be  so  represented  and 
proceeds  to  step  (a)  below.  However,  if  he  is  not  using 
the  beta,  he  enters  the  number  "0"  and  proceeds  to  step 

(3). 

• 

{a)  The  user  specifies  whether  he  will  enter  the 
four-tuples  L,  M,  H,  U  defining  each  beta  variable  direct¬ 
ly  from  the  terminal  or  from  a  data  file  stored  in  the 
computer  by  entering  the  number  "1"  for  data  file  input  or 
the  number  "2"  for  terminal  input. 

(a)  If  he  chooses  the  terminal  input,  his  next 
step  is  to  enter  one  four-tuple  L,  M,  H,  U  for  each  beta 
variable.  If  he  chooses  the  data  file  input,  he  merely 
enters  the  name  of  the  data  file. 

(3)  If  he  is  using  distributions  other  than  the  beta 
to  represent  some  or  all  of  his  component  distributions, 
the  user  enters  the  number  of  distributions  to  be  so  rep¬ 
resented,  followed  by  the  L,  M,  H,  mean,  variance,  and  the 
discrete  cdf  (see  page  C-l)  of  the  distributions  he  has 
chosen.  If  he  is  not  using  other  distributions,  he  enters 
the  number  "0P. 

(4)  If  there  is  a  throughput  (constant) ,  the  user  en¬ 
ters  it  now. 

(5)  He  then  enters  the  number  of  observations  (sample 
size)  he  desires,  followed  by  a  seed  for  the  random  number 
generator. 

(6)  The  computer  prints  the  user's  inputs,  followed 
by  the  output. 

(7)  The  computer  then  queries  the  user  if  he  desires 
another  confidence  interval.  The  user  responds  with  the 
confidence  interval  he  desires,  or  types  the  number  "0"  if 
he  desires  none. 

(8)  The  computer  asks  the  user  if  he  desires  the 
probability  that  some  value  within  the  range  of  the  dis¬ 
tribution  will  be  exceeded.  The  user  responds  with  that 
number,  or  the  number  "0"  if  he  desires  none. 

(9)  The  computer  inquires  to  see  if  the  user  desires 
additional  observations.  The  user  responds  with  the  num- 


ber  of  additional  observations  he  desires  (the  number  "0" 
if  none) .  If  he  desires  additional  observations,  the  com¬ 
puter  repeats  steps  (6)  -  (9) . 

(10)  The  computer  asks  the  user  if  he  desires  to  use 
another  seed  for  the  random  number  generator.  The  user 
responds  with  the  seed  if  he  does,  or  with  the  number  "0" 
if  he  does  not.  If  he  enters  another  seed,  the  computer 
repeats  steps  (6)  -  (10). 
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Program  Listinc 


DOUBLE  PRECISION  A  I  SO  I .H ( 50 ) . XL. XU. Y. ALPHA .UETA 
REAL  L(50l  .H(50>  .MOOE  ISO)  .MM 1 5?))  .U(50)  »C  (56  >  .MEAN (50)  .VAR (50) 
REAL  TITLE (12)  .  VALlNTUO.r.O)  .WI51 ) .PI  I 15) »P2( ID) .P3(15> .P4C15) 
REAL  Cl  US)  .C2tl5l  .C3U5I.C4U5I 

real  MOOE 1 .  ME  AMI .  MOOE  2  .  ME  AN2 .  ME  AN3  >  M00E4  >  ME  AN4  •  MSUM  .LI  .L2.L3.L* 
REAL  WANE  14) 

INTEGER  F1USI.F2U5I.F4U5) 

DATA  NAME/* !EQU* * 'ATE  *.'3  •/ 

11=6  *  • 

JJ=5 
KK=  15 
% 

%  INPUT 


; 


WRITEUI. 5) 

5  FORMA T(//6H  TITLE) 

REAO ( JJ»  7 )  TITLE 
7  F0RMATU5A4) 

WRITE (II .10) 

10  FORMATUH  . 28HNUMBER  OF  BETA  DISTRIBUTIONS) 

REAO ( JJ.  * )  Ml 

IF(N1 .FO.n)  GO  TO  60 
WRITEUI.11) 

11  FORMATUH  .3-’HDATA  FILEU)  OR  TERMINAL  (2)  INPUT) 

REAO ( JJ. * )  INPUT 

IF (INPUT. EO. 2)  GO  TO  16 
HP'.  TEUJ.13) 

13  FORMATUH  .17HNAME  OF  DATA  FILE) 

READ ( JJ. 1 4 )  NAME (4 ) 

14  FORMAT(AU) 

CALL  OUE T (NAME .4 ) 

00  15  1=1. Ml 

REAO  ( 3.  *  )  LU).MOOErI),H(I).UII> 

15  CONTINUE 
GO  TO  60 

j t.  COt* ■  .(Vue 
WRITEUI  .30) 

■30  FORMATUH  .37HL0WEST. MODE. HIGHEST. UNCERTAINTY  COEFF) 

00  40  1=1. N1 
WRITE (1 1. 50) 

50  FORMATUH  4H0ATA) 

READ! JJ.  *  )  L(I)  .MOOEU).M(I).U(I) 

40  CONTINUE 
60  CONTINUE 
WRITEUI. 70) 

70  FORVATUH  .2PHNUM8ER  OF  OTHER  DISTRIBUTIONS) 

REAO(JJ.»)  N2 
IF(N2.EQ.O)  GO  TO  78 
N0>C=N1*1 
NTW0=N1*N2 
WRITEUI. 75) 

75  FORMATUH  .46HL0WEST. MOOE. HIGHEST. MEAN. VARIANCE. DISCRETE  CDF) 
00  76  1=M0CE.NTw0 
WRITEUI. SOI 

REAO(  JJ.*I  L(I )  .MOOE  U )  .HU  )  «MEAN(  I )  .  VAR ( I )  .  (VALINTI J. II  . J=1 . 101 
7 n  CONTINUE 
GO  TO  79 

78  CONTINUE 
NTW0=N1 

79  WRITCUI.72I 

72  FORMATUH  ,  10HTHROUGHPUT) 

REAO ( JJ. • I  TPJT 
WRITE (II. 77) 

77  FORMATUH  .2.>HNUM0ER  OF  OBSERVATIONS) 

REAO ( JJ. * )  M 
WRITEUI.flO) 

SO  FORMATUH  . Z8HSEED  RANDOM  NUMBER  GENERATOR) 

REAO ( JJ. *  1  KSO 
KS02=KS0 

IF(Nl.EO.O)  GO  TO  105 

% 

%  COMPUTE  A  AND  B 
* 

DO  85  lsl.Nl 

N=0  -  - 


I 


NMim  <mooem  )-u<  it  i/tHtn-Min 

V=<  .2B86750fl»U(l )  »••». 

bcii=i.o 

86  ACIt=HII>*NM(Il/(l.-NMtIll 

VUETA=|  (All  Ml.  )»  (HI  l  )M.  )  l/(  (  <A(I)4B(l)42.)*»2.)*(A(I)4BtI)43.l  1 
Z=V-VM£TA 

IFIN.CO.lt  GO  TO  87 
IF(2I  3.86,1 

87  IFIZI  4,85,05 

3  H(Il=0(Il4l. 

GO  TO  86  • 

1  B(lt=B(II-l. 

N=1 

GO  TO  86 

4  Dl 11=0(114.05 
GO  TO  86 

85  CONTINUE 

K 

%  COMPUTE  DISCRETE  COFS 

« 

XL=0.D0 
DO  90  J=1.N1 
XU=1 .00 
ALPHArAt J1 
BETA  =8  ( J  t 

CALL  0QG32  <  XL, XU, T. ALPHA .BETA! 

ccjt=i./r 

xu=o.oo 

VALIMTMO,  J)=l. 

DO  100  1=1.9 
XU=XU4.;00 

call  Ol,G32(  XL.  XU.  V.ALPMA.BETAl 
VAL1 NT ( I «  J 1 =C ( J ) *  Y 
100  CONTINUE 
90  CONTINUE 

% 

*  COMPUTE  CLASS  INTERVALS 

» 

105  CONTINUE 
XLOksTPUT 

xhigh=tput 
00  110  1=1 ,NTwO 
XL0.V=XL0W4L(II 
XHIGH=XrlIGH4H(I  1 

no  continue 

RANGE=XHIGH~XL0W 
W10TM=RANGE/KK 
«r(ll=KLOW 
*(KK41 )=XHIGH 
00  120  I =2, KK 
«r(I)  =  «(I-lt4WlOTN 
120  CONTINUE 

* 

*  GENERATE  HISTOGRAM 

* 

140  CONTINUE 
TOTAU1=0. 

TOTA,.£  =  o. 

TOTAL4=0. 

SUM] =0. 

SUM2=0. 

SUM4 =0. 

SUMS01=0. 

SUMS02=0. 

SUMSO4=0. 

mstar=m 
00  150  1=1. KK 
F  1  ( 1 1  =0 
F2(It=0 
F4 ( I t  =0 
150  CONTtMUE 
160  CONTINUE 
00  180  K=1.MSTAR 
00  |40  J=I.NT»0 
8AN0=U0RNRT IKSOt 


’  *'  it- — -•  -V 


IF1J.GT.1IG0  in  I 95 
T)0  70  9  II..71.IITM0 

00  70S  Irl.io 
MM=I 

IF  (RAnO.LC . VAL1NT 1 1 *LL 1 1  GO  TO  710 
70S  CONTINUE 

710  lFlKM.EO.il  GO  TO  715 

T2=<RAN0-VALIIITIM  u-l.LLl l/I VAL INT  IM.M.LL)  -VALINT  (MM-1  ,LL)  1 
S=IM --1 1/10. 

SMALL 2=.1*T?*S  •  • 

GO  TO  770 

715  SMALL?  —  .  1  *  1RAM0/VALINT1I  >LL1  1 

770  T0TAL?7T0TAL?«SMALl.?*  IHILLl-LlLLl  14L1LL1 

700  CONTINUE 

T0TAL?=T0r.\L?+T7UT 

SUM?7SUM?+ TOTAL? 

SUMS0?7S1J"SQ?-*  TOTAL?*  >7. 

COST?:  (  TOTAL?-  XLO//1 /RANGE 

J?7C05T?4«'«i  . 

IF ( J2.EQ.K -*1 . 1  o?7KK 

•  F?(J2)=F?(J21*1. 

TOTAL270 . 

1<J5  00  2TS  171,10 
M  1=1 

IF(RAN0.LE.VALINT(I,J1 1  GO  TO  ?? 7 
225  CONTINUE 

2?7  IFIM  i.EO.ll  GO  TO  ?31 

2?9  T17IRAN0-VAHNTIM  -1. Jl  1/1  VALINTIMI.JI-VALINTIM  1-1,01  I 
S7(M  -.-li/io. 

SMAL- 1 =. 1*T 1*S 
GO  TO  731 

231  SMALL  1 7. 1* (RAN0/VAL1MTI 1 . J) 1 

?33  TOTAL  1  =  TOT AL 1  *  SMALL  1*  IHIJl-LI  J|  ML1  J1 

T0TAL47T0TAL4 *RAN0* IH1 JI-LIJ1 ) +L I Jl 

190  CONTINUE 

TOTAL1=TOTAL 1 *T?UT 

TO  TALA  ATGTAi.-  .  V—uT 

sum  j  =su**i  4  iotali 
SUM47SUMU*  TOTAL** 

SUMSOIASU'ISIII*  T0TAL1*»2. 

SlJMSnuASU-’SOu*  TOTALA*  *7. 

COST17  1TO.AL1-XLO.T1/RANGE 
J1=C0ST1*K<*1. 

IF  CJ1  .EQ.K/, *1  1  J17KK 
FI 1J1 I7F1 1 Jl )*] . 

COST** 7  (  TOTAL**  -XL0.T1  /RANGE 
J**7C0ST**4K-4*1  , 

IF  I J4.E0.KK* 1 1  J4  7KK 
4IJ4 )7F4 I J4) 4 1 . 

TOTAL  1=0. 

TOTAL470. 

100  CONTINUE 
% 

»  COMPUTE  STATISTICS 

* 

1C=90 

00  235  1=1 ,K< 

Pi  1 1 1 7FL0-TT IF1  1 1 )  1/M 
P21I I7FL0AT1F71 Ill/M 
P«l(  I  I7FL0AT1F4  I  I  1  1/M 
735  CONTINUE 

CALL  STATtSUMl . M.KK.Pl, W* SUM SO l , MEAN1 . MODE  1 r VAR 1.ST01 1 
CALL  C1IMEAN1, IC.P1 .KK.W.wIOTM.Ll.Ull 
CALL  COFIPl.MEANl.KK.W.Cl.PXll 

CALL  STAT I SUM2 «M .K*(  .P2  •  Wt SU-MSO?*  MEAN2  * M00E2 .  VAR?,  ST02 1 
CALL.  Cl  IMEAN2*  IC.P3.KK.Wi  ..'inTH<L7,  U?  1 
CALI.  C0F1P7.KCAN7.K-..W.C7.PX7) 

CALL  STAT  1  SU>14  , M. KK  ,Pu ,  .< •  SIMS04  . ME AN« .  MODEM  ,  VAR4 , STD4 1 
CALL  Cl IMEAN4 , IC.P4 ,KK.W, WI0TH.L4 ,U4 | 

CALL  C0F1Pm,MEAN4,kk,w,C4,PX«1 
IFIN1.E0.01  GO  TO  75? 

00  2S0  1  =  1*111 

ME AN I Il  =  llA(ll*MIIl*n(ll *L III*HIIl*Lflll/IAIll 4H III ♦?# 1  I 

VARI117|IAI11«l.)*tHlll*l.l*UH(l|-L(I)l«*?.)l/mAllMUIIl*2.  1  * • 2. > * t A ( I >  *0 1 1 1 *3. )  I 
250  CONTINUE 
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4*-  ,  ■  t-NT  ; 


A"  '»'%*• 


o 


226 

227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 

241 

242 

243 

244 

245 

246 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 
261 
Cut 

263 

264 

265 

266 

267 

268 

269 

270 

271 

272 

273 

274 

275 

276 

277 

278 

279 

280 
281 
282 
283 
2B4 

285 

286 

287 

288 

289 

290 

291 

292 

293 

294 

295 

296 

297 

298 

299 

300 


252  CONTINUE 

MSUM=TPUT 

VSUM=0. 

DO  25‘.  I=1.NTW0 
MSUM=M5UM.M£AN  I 1  ) 

VSUMrVSUM.VAR! 1 ) 

255  CONTINUE 

* 

MEAN3=(XHIGH+XL0WI/2.  .  . 

VAR3=<1  ./12.  1  »RAMGE»*2.  * 

ST03=S0RT( VAR3I 
UNI  =  1./Ki< 

00  290  1=1. KK 
P3T  I*=UNI 
290  CONTINUE 

CALL  Cl (MEAN3. IC .P3.KK.W. wlOTH.L3.U3* 

CALL  C0F(P3.MEAN3.K<.W.C3.PX31 

* 

*  PRINT  INPUT 

« 

WRITE  1 1 1 .295* 

295  FORMAT (1H  * 

WRITE!  11.30  1* 

300  FORMAT!/ I 
WRITEI1I.3101 
310  FORMAT!//) 

WRITE (11.320) 

320  FORMAT!////* 

wRI TE ( 1 1 .330 )  TITLE 

330  FORMAT I 1H  .55X.15Att) 

WRITEI1I.320* 

WRITE (I  I. 335* 

335  FORMATdH  . 55X.11HI  N  P  U  T  S> 

WRITEIII.320I 
WRITE! 1 1 .390)  M 

390  FORMATdH  .  15H0HSERVAT10NS  =  .151 
WR Z  7C 1 1 1  f  •*0..'  )  r 

400  FORMATdH  .7HSEE0  =  »5X.  18) 

WRITE (1 1. 310) 

WRITE (II .3401 

340  FORMATdH  .9HRETA  V8LE  .  7  X  ,  6HL0WEST  .  6X.  4HMOOE  .5/.  7HHIGHEST  •  &■ 
5X.7HU  COEFF  »  5X»  5HNV00E  .5X  .  5HALPHA .  6X ,  4HBETA .  7X»  4HMEAN.  5X .  BHVARI  A:*CE> 


WRI TE ( I i . 295) 

1FIN1.EO.O)  GO  TO  351 

WRITE (I  I .350*  ( I .LI  I ) • MOOE I  I ) .HI  I  )  .U(  I  I  .NMI 1 J , A ( I ) .8  I  I ) .MEAN! I ). VAR  II )  . 
350  FORMATdH  .  3X.  I2.AX.F9.]  .  1  X.F9.  1 . 3X.F9. 1 . 6X. F4 .2 . 7X »F4 , 2 .* 
6X.F5.1.5X.F5.1.2X.F9.1 , IX. FI  2. 1 > 


1FIN2.EO.0I  GO  TO  352 
WRITE(II.30ni 

351  CONTINUE 
WRITE! 1 1.353) 

353  FORMATdH  .10HOTHER  VBLE* 

WRITE! II . 354 )  (  I  .L(  1  ) .MODE  I  I  * .H( 1 ) .MEAN!  I  I . VAR ( I* . 1 SNONE.NTwO) 

354  FORMAT  (li!  .3X.  I2.8X.F9. 1 . 1  X.F9 . 1 . 3X.F9,  i  ,miX.F9. 1 ,  1  X.F12. 1 1 

352  CONTINUE 
WP*T£!M  »2°5> 

WRITE (I I. 355*  TPUT. TPUT.TPUT 

355  FORMATdH  . lOHTHROUGHP'.lT. 3X.F9.1.13X.F9.1.44X.F9.1) 


WRITE!! I .3001 

WRITE II 1 .360)  XLO'I.  XHI5H.MSUM.  VSUM 

360  FORMAT! 1H  .8h*.SUMS»»»5X.F9.I . 13X.F9. 1 . 44X.F9. 1 .3X.F1 0. 1 ) 
IFIN2.EO.O)  GO  TO  410 
WRI TE 1 1 1 .300 ) 

WRITE (II .320) 

WRI TE ( 1 1 » 370 ) 

370  FORMATdH  . IOHOTMEr  VRLE . 50 X . 3HCOF * 

WRITE (II. 300  I 

WRITE (I I.3«0)  (I • ( VALINTI J. I  * . J=1 . 1  0) . I=NONC.NTWO) 

380  FORMATdH  <3X.  I2.20X.  1 0  (F6.4 ,3X1 ) 

WRITE!! I .320* 

* 

*  PRINT  OUTPUT 

% 

410  CONTINUE 
WR1TCIII.320I 


1  =  1.  Ill ) 
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.5SX.1IH0  U  T  P  U  T> 


SOI 

302 

303 
30* 

305 

306 

307 
30A 

309 

310 

311 

312 

313 

314 

315 

316 

317 

318 

319 

320 

321 

322 

323 
32* 

325 

326 

327 

328 

329 

330 

331 

332 

333 
33* 

335 

336 


f 


o 


5ENCC  INTERVAL?) 


I 

f 


WRITEM  l.«O0) 

*90  FORMAT  1 1H  »55x» 1 1HO  U  T  P  U  T) 

WRITE  II 1.320) 

WR  I  7£  1 1 1 #500 1 

500  FORMATMH  ,  AH  INTERVAL  » 1  OX ,  5MRANGE  »  6X, 22H INDEPENDENT  nETA/OT)CH.3X.» 

20HOEPLIIUENT  IIETA/OTHER.3X.13HTOTAL  UNIFORM. SX»  19HIN0EPF.N0ENT  UNIFORM) 

WRITE  1 1 1. 295) 

WRITE  Cl  I. 505) 

505  FORMAT  UK  .36X,3HPnF,ux,3HCDF.12x.3HP0F.4X.3HC0F.10X.* 
3HP0F.4X.3HCDF.12X.3MPDF.4X.3HC0F)  < 

WRITE) II .300) 

WRITE  1 1 1.510)  f I . Wf I ) .  wt  1*1 )  .PI  1 1 )  ,C1  II )  .P2CI ) .C2CI ).P3CI ) #C3C 1 1 .P* C I ) »C* C I ) » 1  =  1 •KW) 
510  FORMATMH  .2X. I2.4X.F9.1 ,3H  -  ,F9.  1 .5X.F6.4, 1  X.F6.4 .9X.F6.4  .  1  X.F6.U, 7X.* 
F6.4,2X.F6.4,RX,F6.4,2X.F6.4) 

WRITE  til. 310) 

WRITECI  I  .520)  MEAN1  .**E  AH2  .  MEACI3  .MEACI4 

520  FORMATMH  .Ah*  *MEAt|..,p?x,F9. 1  .  13X.F9.1  .  I2X.F9. 1 . 1 3X.F9.  1 ) 

WR  I  TE  1 1 1 . 530  )  MODE  1  .  MOI  )E2  .  MODE  4 

530  FOOMATMH  .«•)». MODE*  ■  .27X.F9.  1  .  13X.F9.1 .34X.F9.I) 
wRITEC I  I .540  I  VAR  t . VAR2.  VFR3, VARU 

54  0  FORMAT  MH  . 1  PH*- VARIANCE* ‘ .20X. FI  2.1 .10X.F12.1.9X.F12.1.10X.F12.1) 

WRITECI I  .545)  STOl  .ST02.STD3.ST0U 

545  FORMAT  C 1 H  .17H**STD  DE VI A TION* *. 1 7X . FI  0 . 1 . 12X.F1 0. 1 . 1 1 X.FlO. 1 . 12X .F 1 0. 1 ) 

WRITE  1 1 (.295) 

620  CONTINUE 

WRITECI?. 550)  IC.L1.E2.L3.L4.U1.II2.U3.U* 

550  FORMATMH  . 2H *  * ,12.2 3HV  CONFIDENCE  INTERVAL** .RX.F9.1 . 1H. . I2X.FR. 1.1H. .V 
11X.F9.1. iH. .12X.F9.1 .1H./36X.F9.1.13X.F9.1.12X.F9.1.13X.F9.1) 

1FC7C.NE.90)  GO  TO  571) 

WRITECII.295) 

WRITECI), 560)  oxl.PX2.PX3.PX4 

560  FORMATC1H  .20H*«PROR  EXCEED  ME AN* * .22X.F3 .2 , 19X.F3. 2 . 1 BX.F3.2 . 19X.F3. 2) 

WRITECI I .300) 

570  CONTINUE 
WRITE)! J. 310) 
w«ITE  MI. 600) 

60-)  FORMAT:::!  . 2S:'.ANCT) :ED  C: 

RE  AO  C  J-  J,  *  )  1C 
IFCIC.EO.O)  GO  TO  610 
CALL  Cl  cmeani  .  IC.P1  .K-f  ,w,  width.  LI  ,U1 ) 

CALL  Cl  IMEAN2.  IC  ,P2.K^  ,  w,  IDTH.  L2  •  U? ) 

CALL  Cl  CMEAN3.  IC  . P3 •  Kx ,  >•;.  .vIOTH.L3.U3) 

CALL  C  1 -CMEAN4  «  IC  ,  P4  ,  K  \  ,  vt,  VIOTH.L4.UU) 

GO  TO  t, 20 
610  CONTINUE 
WRITECII.310) 

WRITEtlI.622) 

622  FORMATMH  .23HPROB  EXCEED  SOME  VALUE?) 

REAOCJ.i,*)  ? 

1FCZ.E0.0. )  GO  TO  623 
CALL  COF  CPI .  Z.Kx.  W.C1  .PXI  ) 

CALL  C0FtP2,2.KK,w.C2.PX2) 

CALL  CC)FIP3.Z.KK,W,C3,PX3) 

CALL  CDF  CPU , Z.WK , W.CU ,PX4 ) 

WRITE ( 1 1 .624)  Z.PX1.PX2.PX3.PX4 

624  FORMAT 1 1H  ,13H»*PROf)  EXCEED ,F 1 0 . 1 . 2H* * ,  17X.F3.2, 19X.F3.2, 1BX.F3.2 . 19X. F3. 2) 

GO  TO  6)0 

623  CONTINUE 
WRITECII.310) 

WRITE  1 1 1 .630) 

630  FORMATMH  .24HA0l.'ITI0NAL  OBSERVATIONS?) 

READ ( Jj. *  1  MSTAR 
IFCMSTAR.EO.O)  GO  TO  640 
M=M* MSTAR 
GO  TO  160 
640  CONTINUE 
WRITECII.310) 

WRITE  II 1.650) 

650  FORMATMH  ,  1 3HANOTHER  SEED?) 

REAOCJJ.O  KSO 
IFCKSU.EO.O)  STOP 
KSD2=KSO 
WRITCCII.77) 

REAOIJJ.*)  M 
00  TO  1*0 


C-12 


^  Jr~f  ZT-.-r 


© 


1  SUBROUTINE  Cl  (MEAN,  1C, P,KK,W.  WIDTH, LBU.tJttO) 

2  REAL  w ( 1 ) , Mt AN, P <11 , LwTl , LPROB, LBD ,LWT2 

3  DO  2 SO  1=2, KK 

4  J=J 

5  1F(W< I )  .GE.MEAN)  GO  TO  260 

6  2S0  CONTINUE 

7  260  MM=0 

8  CON=.5*(IC/100.» 

9  LWT1=(MEAn-w( J-l ) >/( w( J)-w(U-l) > 

10  270  LPR0B=LWT1»P(J-1) 

11  K=1 

12  IP ( LPROO, GE • CON)  GO  TO  290 

13  '  J1=J-1 

14  DO  280  1=2, J1 

15  K=I 

16  LPROB=LPROB+Pt J-II 

17  IF (LPROB. GE. CONI  GO  TO  290 

18  200  CONTINUE 

19  CON=2.«CON-LPROB 

20  LBD=V;(1I 

21  GO  TO  300 

22  290  LWT2=(LPROB-CON)/P< J-K) 

23  LBD=(.'<J-K)+LwT2*wI0TH 

24  IFIMC.EO.I)  GO  TO  3  0 

25  S 

26  300  RwTl=l.-LWTl 

27  RPROO=R*T1*P< J-l) 

28  LL=1 

29  IFtRPROB.GE.CON)  GO  TO  320 

30  J2=KK-J»2 

31  DO  310  1=2, J2 

32  LL=I 

33  RPROB=RPROB+P<J+I-2l 

34  IF <RPROB.GE.COM)  GO  TO  320 

35  310  CONTINUE 

36  CON=2.*CON-RPROB 

37  U0U=w (KK+ 1 I 

38  KM=1 

39  GO  TO  270 

40  320  RwT2=(RPf<0B-C0N)/P<  J4LL-2I 

41  UBO=«r<J+LL-lI-RwT2»WIDTH 

42  330  RETURN 

1  SUBROUTINE  CDF|P,z,KK,*rC,PROBXJ 

2  REAL  P(1I«V<1>*C(1I 

3  C(1I=P(1) 

4  DO  10  1=2,  Kl. 

5  C(l>=i-<1-1)+Pl  J) 

6  .  10  CONTINUE 

7  N=KK41 

8  00  20  1=2* N 

9  J=1 

10  IFIWdl.GT.Zl  GO  TO  30 

11  20  CONTINUE 

12  30  A=<z-wtj-m/<w<j»-wij-i>) 

13  PROBX=(l.-A)*P<U-l) 

14  DO  40  I = J, KK 

15  PROBX=PROUX+P ( 1 ) 

16  40  CONTINUE 

17  RETURN 
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1  SUBROUTINE  00G321XL.XU. Y. ALPHA . BETA) 

2  DOUBLE  PRECISION  XL. XU» Y. A. B.C «FCT» ALPHA. BETA 

3  FCT(X)=X»-.ALPHA*(I.-X)*»B£TA 

4  A=.5D0*(XL*XU) 

5  B=XU-XL 

6  C= • 49863 1930924 740 7800 *B 

7  T=(. 3509305004 73504630-2 >*(FCT(A*C>«FCT( A-C >  > 

6  C=. 4928057b  7/26341700*8 

9  Y=Y+ (.  61371973654528350-2)  *t  PC  T  l  A+CUFCT(A-C>  ) 

10  C=. 4823811277937532200*6 

n  Y=Y* ( .126960  326546310300-1 1  * ( FCT I A*C ) *FCT ( A-C ) I 

12  .  C=.46745303796Bco934D0*8 

13  Y=Y*(  .  171 36931 45651 0  71  70-1)  *  C  FCT  ( A*C  ) +FCT  t  A-C)  ) 

14  C=. 4481605778630260600*8 

15  Y=Y-M  .214179490111 ' 33400-1 1 *( FCT ( A*C ) +FCT C A-C > ) 

16  C=. 4246838 068 6c 2849  D0*8 

17  Y-  Y*  ( .  25u9  0296311  3'  OP..O-l)*tFCT(A+C)*FCTIA-C) ) 

16  C-. 3972 41 897983971 2OD0*b 

19  Y=Y* ( .2934  204  6739267 7  740-11* l F CT ( A+C ) +FCT C A-C I  1 

20  Cf. 3660910693701448400*8 

21  Y=Y+(. 32911  1113861  -.09230—1  >*(FCT(A*C)*FCT  1  A— C 1 1 

22  C-. 3315221334661076000*8 

23  Y- Y* ( . 361 72897054 -242530-1 )*(FCT(A*C)+FCT l A— C 1 ) 

24  C=. 2936576786203611600*8 

25  Y=Y*(. 3909694  78935351 530-1 1  * ( FC T ( A*C i +FCT 1 A-C) ) 

26  C=. 2534. 94-54466  114/000*8 

27  Y=Y* (. 616559621 134733780-1 )* (FCT I A*C)+FCTt A-C) ) 

C-.  21 067563  *06621 76700*5? 

29  Y=Y*( .4362604650220 19060-1 )• (FCT (A*C)*FCT( A-C > ) 

30  C=. 1659343011410638200*8 

31  T=Y*  (.455869392478-'19420-l)*(FCT(A+C)*FCT(A-C») 

32  C=.l 19643681 1260685400*8 

33  Y=Y* ( . 469221 9-5404022830-1 ) * (FCT (A+C ) *FCT ( A-C) ) 

34  C=. 72235980791393250-1*8 

35  Y=Y*(  ,47819360  '396374300-1)*(FCT(A*C)*FCT(A-C1  .• 

36  C=. 241533328438691580-1*8 

37  Y=B* ( Y* ( .482700442573639000-1 ) * (FCT (A+C ) +FCT (A-C) ) ) 

38  RETURN 

39  END 


S  SUBROUTINE  STAT ( SUM. M.KK*  P. B. SUMSQ. ME AN. MODE. VAR. STD) 

2  REAL  MEAN.MOOE.Plll.YKU 

3  KEAN=SUM/M 

4  MOOE=0, 

5  DO  10  1=1. KK 

6  IFtP(I).GT.MOUE)  MODE=P(I> 

7  1FIP(I>. Ell. NODE)  I  MODE  =  I 

0  .  10  CONTINUE 

9  MOOE= •  5*  (  w  ( I  MODE  )  +  W  ( 1  MODE*  1 )  ) 

10  VAR=(SUMSQ-M*N£AN*»2. >/(M-l> 

11  STO=SQRT ( VAR) 

12  RETURN 
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Optimal  Sample  Size 

In  most  sampling  procedures  the  larger  the  sample  the 
closer  the  sample  distribution  approximates  the  true  dis¬ 
tribution.  But  larger  samples  are  more  expensive  to  gen¬ 
erate  than  smaller  ones.  The  simple  experiments  described 
below  represent  an  attempt  to  determine  an  optimal  sample 
size  for  Program  SPET  —  the  smallest  sample  size  that 
will  ensure  "reasonable"  accuracy  in  the  sampling  proce¬ 
dure  . 


A  statistic  called  the  K  statistic  in  this  study  is 
used  in  the  search  for  the  optimal  sample  size.  It  is  de¬ 
fined  as 


K  =  E  (Oj  -  ei)2 
i=l  1  1 


where  N  =  the  number  of  class  intervals 


0£  =  the  number  of  observations  occurring  in  the 
ith  interval  *  M. 


M  *»  -he  total  number  of  observations 

e^  =  the  number  of  observations  expected  in  the 
ith  interval  (given  that  the  process  gener¬ 
ating  the  observations  is  following  a  par¬ 
ticular  statistical  distribution)  *  M. 


The  K  statistic  is  the  sum  of  squared  deviations  of  the 
observed  probabilities  from  the  expected  probabilities  of 
each  class  interval.  As  the  size  of  a  randomly  drawn  sam¬ 
ple  is  increased,  the  K  statistic  decreases  in  value  until 
lim  K  =  0. 

M  ® 


The  first  experiment  consists  of  using  five  randomly 
selected  seeds  with  the  uniform  random  number  generator 
used  by  Program  SPET  to  generate  five  sequences  of  K  sta¬ 
tistics.  Each  sequence  contains  a  K  statistic  for  sample 
sizes  500,  1000,  2000,  ...,  10000.  These  K  statistics  are 
plotted  in  Figure  C-l  on  page  C-16.  Note  how  the  se¬ 
quences  converge  at  sample  size  6000.  It  appears  that 
this  might  be  the  optimal  sample  size.  Can  one  expect  a 
sample  size  of  6000  to  ensure  "reasonable"  accuracy  in  the 
sampling  procedure? 
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FIVE  SEQUENCES  OF  K  STATISTICS  AS  A  FUNCTION  OF  SAMPLE  SIZE 


1 


The  accuracy  of  the  Monte  Carlo  sampling  procedure 
used  in  Program  SPET,  for  purposes  of  this  inquiry,  is 
measured  in  terms  of  the  percent  deviation  of  certain  sta¬ 
tistics  from  their  true  values.  The  second  experiment  is 
an  attempt  to  measure  the  accuracy  of  Program  SPET  for 
various  sample  sizes.  It  consists  of  using  the  same  five 
seeds  selected  in  the  first  experiment  to  draw  five  se¬ 
quences  of  samples  from  a  uniform  distribution.  Each  se¬ 
quence  contains  samples  of  size  500,  1000,  2000,  3000, 
6000,  and  9000.  The  mean  and  the  lower  and  upper  confi¬ 
dence  limits  of  the  90%  confidence  interval  are  noted  from 
the  output  of  Program  SPET  and  the  percent  deviation  of 
these  statistics  from  their  population  values  is  computed. 
Then  the  maximum  of  the  absolute  value  of  the  deviations 
is  selected  for  each  statistic  in  every  sample  size  and 
plotted  in  Figure  C-2  on  page  C-18.  Note  that  the  rate  of 
decrease  in  the  error  (maximum  percent  deviation)  of  these 
statistics  is  rapid  in  the  range  of  the  sample  sizes  500 
to  2000,  slowing  somewhat  after  sample  size  2000. 

Consider  the  error  in  sample  size  2000.  Would  the 
expectation  of  a  deviation  of  at  most  .21  percent  in  the 
mean  and  2.44  and  1.81  percent  in  the  lower  and  upper  lim¬ 
its  of  the  confidence  interval  respectively,  be  "reason¬ 
able?"  The  authors  would  answer  affirmatively.  Reason¬ 
ableness  is  subjective.  It  is  felt  that  the  accuracy  of 
sample  size  6000  is  not  enough  better  than  that  of  sample 
size  2000  to  warrant  incurring  the  increased  cost  of  gen¬ 
erating  an  additional  4000  observations. 

Much  greater  confidence  could  be  placed  in  these  ten¬ 
tative  observations  if,  instead  of  five  sequences,  30,  40 
or  more  sequences  had  been  generated.!./  But  even  the  re¬ 
sults  of  the  five  sequences  permit  a  more  confident  choice 
of  sample  size  than  no  experimentation  at  all. 


£/  Along  with  a.  gh.ia.tiK  numbeK  AequenceA  one  might 
have  Kepeated  expeximent  two  u&ing  one  ok  two  xepxe- 
Aentative  beta  diAtxibution A  in  addition  to  the  uni- 
£oxm  diAtKibution  uAed  above. 


C 


SAMPLE  SIZE 


MAXIMUM  PERCENT  DEVIATION  (ERROR)  OF  CERTAIN  STATISTICS 
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